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The breakthrough capacity of optoacoustics for three-dimensional visualization of dynamic events in real 
time has been recently showcased. Yet, efficient spectral unmixing for functional imaging of entire 
volumetric regions is significantly challenged by motion artifacts in concurrent acquisitions at multiple 
wavelengths. Here, we introduce a method for simultaneous acquisition of multispectral volumetric datasets 
by introducing a microsecond-level delay between excitation laser pulses at different wavelengths. Robust 
performance is demonstrated by real-time volumetric visualization of functional blood parametrers in 
human vasculature with a handheld matrix array optoacoustic probe. This approach can avert image 
artifacts imposed by velocities greater than 2 m/s, thus, does not only facilitate imaging influenced by 
respiratory, cardiac or other intrinsic fast movements in living tissues, but can achieve artifact-free imaging 
in the presence of more significant motion, e.g. abrupt displacements during handheld-mode operation in a 
clinical environment. 



First optoacoustic images from living biological tissues were demonstrated over a decade ago 1 ' 2 . The unique 
advantages of this hybrid imaging technology drive continuous technical advancements as well as its use in a 
number of novel biomedical applications. In optoacoustics, sound waves are induced by absorption of time- 
variant light radiation 3 , which allows for high-resolution imaging of rich optical contrast beyond the depth- 
penetration limits imposed by light scattering. Thus, otherwise unresolvable or undetectable anatomy and 
function can now be visualized at millimeter to centimeter scale depths in living organisms 4 ' 5 . 

Translation into clinical practice occurs to be the natural next step so that the unique strengths offered by the 
optoacoustic technology can be fully exploited. Several dedicated systems using endoscopic, hand-held, or 
stationary design approaches have been recently proposed for clinical use, with multiple applications envisioned 
in diagnosis of breast cancer 6,7 and cardiovascular disease 8 , ocular imaging 9 , visualization of arthritis 10 , lymphatic 
system 11 , skin abnormalities 12 , prostate 13 and thyroid glands 14 . Of particular convenience is the use of hand-held 
probes, which permits real-time visualization of dynamic events, and has the potential to speed up clinical 
examinations on a broad range of applications. In fact, these are also the main assets that propelled ultrasono- 
graphy into the most- widely- employed imaging modality in today's clinical practice. On the other hand, the 
optoacoustic tomographic reconstruction problem and physics behind the signal excitation substantially differ 
from ultrasonography. While linear arrays are typically used in ultrasound imaging, curved array geometries 
proved to be a significantly better strategy for collection of optoacoustic responses and further resulted in more 
accurate reconstructions 1516 . Moreover, optoacoustic excitation via laser light offers additional exclusive capa- 
cities, such as intrinsic ability to acquire full tomographic datasets within a single interrogation laser pulse, thus 
achieving real-time three-dimensional imaging capacity 17 " 19 . Most importantly, by acquisition of images at 
multiple optical wavelengths, multispectral optoacoustic tomography (MSOT) can map functional contrast as 
well as the bio -distribution of molecular agents in tissues 20,21 . 

Real-time acquisition of multispectral datasets imposes significant technological challenges when considering 
motion in living subjects, e.g. due to breathing or heartbeat. Even more significant artifacts may be introduced by 
operating optoacoustic imaging devices in a hand-held mode. Several prior attempts to avoid motion-related 
artifacts in multispectral optoacoustic imaging were done using high repetition lasers in optoacoustic endoscopy 22 
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and optical -resolution photoacoustic microscopy 23 , as well as by 
means of Raman shifters in acoustic resolution optoacoustic micro- 
scopy 24 . Yet, these implementations employ a single detector 
approach, which only renders real-time performance in two dimen- 
sions and is overall not suitable for hand-held operation. Most 
recently, a true five-dimensional (5D) optoacoustic tomographic 
approach has been demonstrated by fast laser wavelength tuning 
and real-time acquisition and visualization of three-dimensional 
multispectral data using spherical matrix array detection geometry 25 . 
Yet, the substantial (millisecond level) time delay between the excita- 
tion laser pulses at different wavelengths has severely limited the 
application in hand-held mode. Significant motion artifacts intro- 
duced by hand-held scanning were observed when the velocity of the 
probe exceeded several mm/s. The displacement d of the probe 
between consecutive laser pulses is given by d = v/PRF, where v is 
the scanning velocity and PRF is the pulse repetition frequency of the 
laser. For instance, considering a pulse repetition frequency of 50 Hz 
and a relatively slow hand-held scanning velocity of 5 mm/s, the 
probe moves by 100 um between two consecutive pulses at two 
different wavelengths, which prevents accurate spectral unmixing 
preserving high spatial resolution. Motion effects scale with the num- 
ber of wavelengths used for multispectral data acquisition, further 
affecting the accuracy of real-time multispectral imaging. 

In this work, we explore the feasibility of functional volumetric 
(three-dimensional) optoacoustic imaging in a hand-held mode by 
introducing a microsecond-range (as opposed to millisecond- range) 
delay between laser pulses at different wavelengths. In this way, the 
time lapse required for collection of multispectral data can be mini- 
mized, eliminating potential motion artifacts and enabling func- 
tional multispectral imaging even during the presence of abrupt 
displacements of the imaging probe. 

Results 

Experimental concept. The layout of the experimental system 
employed is depicted in Fig. 1. Acquisition of three-dimensional 
tomographic data in real time was performed by an optoacoustic 
hand-held probe consisting of an ultrasound array with 256 
piezoelectric elements distributed on a spherical surface 26 . Details 
on the array geometry are provided in the methods section. 
Optoacoustic responses in the imaged object are excited by 
nanosecond- duration pulsed beams from two tunable OPO lasers 
(Phocus, Opotek Inc., Carlsbad, CA), which were set at wavelengths 
of 760 nm and 850 nm. The light was guided through a custom- 
made fiber bundle (CeramOptec GmbH, Bonn, Germany) having 
two input arms and a common output arm. The individual fibers 
at the output of the bundle were randomly distributed with respect to 
the two inputs so that both beams created approximately the same 
illumination profile at the tissue surface. The first laser was internally 
triggered whereas the flash lamps and Q-switch of the second laser 
were externally triggered and conveniently delayed as shown in 
Fig. 2a. Two delay generators were used in the synchronization 
scheme. The first delay D x was selected according to D 2 = T-L 
from the Q-switch output of the first laser, where T and L denote 
the pulse interval and the resulting time lapse between the two laser 
pulses, respectively. Since both corresponding time-resolved signals 
were detected within the same acquisition window, the time lapse L 
needs to be sufficiently long to separate the reconstructed images. 
Given a pulse interval of T = 100 ms, a value of 99.983 ms was 
chosen for D x which translates into a time lapse of L = 17 |is 
between the Q-Switch trigger signals of the two lasers. This time 
lapse corresponds to approximately 26 mm of propagation of the 
ultrasound wave within tissue, which is larger than the imaged region 
of interest, thus no crosstalk between two consecutive images 
acquired at the two wavelengths is expected. Since the second laser 
is triggered externally, both Q-Switch and flash lamps needed to be 
synchronized. To ensure optimal signal energy, a delay of C = 180 jas 
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Figure 1 | Lay-out of the experimental set-up. The computer picture was 
developed by the authors. QS-Q-switch, FL-flash lamp. 

had to be set between the respective trigger signals. This was 
maintained by setting the delay of the second delay generator D 2 
to a value of 99.820 ms from the Q-switch input of the second 
laser according to D 2 = T-C. 

Imaging of tissue vasculature in a healthy volunteer. The palm of a 
healthy volunteer was scanned in order to test the in- vivo imaging 
performance of the introduced method. An example of optoacoustic 
signals acquired from two microsecond-delayed laser pulses is 
displayed in Fig. 2b. It can be seen that the optoacoustic responses 
due to pulses at two different wavelengths can be separated on a time 
axis without cross talk. Yet, the short delay of only 17 |is between the 
two pulses ensured that no motion artifacts affect the unmixing of the 
multi- spectral data. The imaging results are showcased in Fig. 3 while 
two additional movie files are available in the online version of the 
manuscript that best illustrate the attained imaging performance in a 
dynamic fashion. Image reconstruction and unmixing was 
performed as described in the methods section. Figs. 3a and 3b 
display examples of maximum intensity projection (MIP) images 
acquired respectively at 760 nm and 850 nm, both color-coded for 
depth. For the same imaged volume, images acquired at 850 nm 
show a different appearance than those at 760 nm, which is mostly 
related to the strong wavelength-dependent absorption by the blood 
chromophores. For instance, the highly oxygenated vessels (arteries) 
are naturally expected to present higher signal intensities for the 
850 nm wavelength as compared to 760 nm, whereas an opposite 
effect is expected for less oxygenated vessels (veins). The distinction 
between vessels with different blood oxygenation levels is best 
perceived in the blood oxygen saturation maps displayed in Fig. 3c. 
One should notice that wavelength-dependent light attenuation may 
lead to uncertainties (errors) in estimating the blood oxygen 
saturation, which also scale with depth. Our analysis, based on the 
procedure detailed in the methods section, reveals that for vessels 
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Figure 2 | (a) Synchronization scheme between the two lasers, 
(b) Example of the typically acquired optoacoustic signal train. 
QS - Q-switch, FL - flash lamp, DAQ - data acquisition system. 

shallower than 6 mm (as shown in Fig. 3), such uncertainty is 
generally expected to be lower than 15% (Fig. 4). However, 
significant errors may be introduced at deeper locations and a 
model accounting for the effects of light attenuation in tissues is 
needed for accurate estimation of the oxygen saturation. The 
dynamic imaging capacity is further showcased in the two movie 
files available in the online version of the journal. In this case, the 
MIPs of the unmixed distribution of oxygenated and deoxygenated 
hemoglobin are displayed in arbitrary units in all three axes for better 
distinction of vessels having different oxygen saturation levels. Movie 
1 demonstrates the ability to resolve distribution of oxygenated and 
deoxygenated hemoglobin in a pulsating artery from the data 
captured by the resting hand-held probe system. In Movie 2 the 
system is shown to attain dynamic tracking of deep human 
vasculature in real time by moving the handheld probe along the 
skin surface. In both cases, no motion artifacts were detected in the 
images corresponding to the different wavelengths, so that fast 
acquisition of multispectral datasets enables generation of motion- 
artifact-free functional images. 

Comparison to the fast wavelength tuning approach. To compare 
the performance of the suggested method to the fast wavelength 
tuning approach, we simulated the latter by using instead the 
multispectral data from two adjacent laser pulses separated by 
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Figure 3 | Optoacoustic images of tissue vasculature in a healthy 
volunteer acquired using the three-dimensional handheld probe and the 
developed microsecond delay approach, (a) Volumetric images at 760 nm 
color-coded for depth, (b) The corresponding images at 850 nm. (c) Blood 
oxygen saturation maps (maximum intensity projections along the depth 
direction) derived by spectral unmixing of the multispectral data. The 
scalebar corresponds to 5 mm. 

approximately 100 ms. The reconstructed MIP representations of 
the volumetric images acquired at 760 nm and 850 nm using both 
methods are compared in Fig. 5. A manual (feature-based) motion 
correction of the images at the two wavelengths has been further 
performed for the simulated fast wavelength tuning approach. 
Motion- related errors were not identified in images acquired using 
the microsecond-delay data (panels a-c). However, motion-related 
errors were prominent in images acquired with the fast wavelength 
tuning approach (panels d-f). The motion between images can be 
best recognized when superimposing the images at 760 nm onto the 
respective images at 850 nm (see the green squares in Figs. 5b and 
5e). Clearly, when the two images are separated in time by 100 ms, 
the structures are shifted by approximately 800 um between adjacent 
pulses at the two different wavelengths, making spectral unmixing an 
impossible task during even a relatively slow scanning speed of 
8 mm/s. Manual correction for motion before performing the 
unmixing step is still not able to rectify all the motion-related 
artifacts in the unmixed images (Fig. 5f). For example, the vessel 
marked 'A' obviously generates similar signal intensities at the two 
wavelengths when the data is acquired using the microsecond- delay 
approach (Figs. 5a and 5b), whereas the same vessel generates a 
stronger signal level at 760 nm versus 850 nm when the acquisi- 
tion at the two wavelengths is separated by —100 ms (Figs. 5d and 
5e). This discrepancy can be attributed to a shift in the light fluence 
and sensitivity maps of the detection array as it moves to a slightly 
different location. As a result, unmixing of motion-affected data 
would give rise to erroneous interpretation of low oxygen 
saturation in the vessel (Fig. 5f versus Fig. 5c). The opposite effect 
occurs for the vessel marked 'B', which appears to generate much 
stronger signal intensity at 850 nm due to motion. As a result, the 
blood oxygen saturation in this particular vessel is overestimated. 

Discussion 

The results demonstrate the feasibility of three-dimensional real- 
time functional imaging free of motion -related artifacts. The pre- 
sented approach was demonstrated to successfully mitigate artifacts 
in the presence of intrinsic tissue motion, such as circulation or 
respiratory movements, but also when fast scanning with a handheld 
probe is performed. For the dual wavelength approach used in the 
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Figure 4 | Uncertainty (error) in approximating the blood oxygen saturation s0 2 , a of a target point located at different depths by assuming 
wavelength independent light fluence 0(r) in Eq. (2). Three different oxygen saturation values of the background tissues s0 2 ,b are considered, namely 
90% (a), 70% (b) and 50% (c). 



current study, even moving the probe at unreasonably high velocity 
of 2 m/s with respect to the imaged object would only result in 
displacement of 34 (im within the 17 |is delay between laser pulses 
at the two wavelengths, which is below the achievable spatial resolu- 
tion of the three dimensional imaging probe. Furthermore, using the 
proposed method, hand-held scanning can be performed with low 
laser repetition rate of several tens of pulses per second allowing for 
high per-pulse energies in the tens and hundreds millijoule range that 
still conform to laser safety standards, thus guaranteeing good signal- 
to-noise performance and deep tissue penetration. For comparison, 
in order to reduce motion artifacts in conventional multispectral 
optoacoustic tomography using a single laser source, the repetition 
rate of the laser must be significantly increased in order to achieve 
fast wavelength tuning and minimize the time lapse between pulses at 
different wavelengths 25 . However, this would not constitute a viable 
alternative to the approach suggested here since increasing the pulse 
repetition rate would readily imply corresponding reduction of the 
per-pulse energies in order to conform to the laser safety standards 
on the average light intensities delivered to the imaged subject 27 . 
Decreasing the per-pulse energy would come at the adverse 
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Figure 5 | Maximum intensity projections of the images acquired with 
two laser pulses separated by 17 /us (a-b) and the corresponding image 
representing the unmixed distributions of oxygenated and deoxygenated 
haemoglobin (c). The corresponding images acquired with two adjacent 
laser pulses separated by approximately 100 ms are shown in (d-f). Part of 
the image data in (a) and (d) is superimposed in green into the areas 
marked with the green squares in (b) and (e). Maximum intensity 
projections along the depth direction are shown. The scalebar corresponds 
to 5 mm. 



compromise of reduced signal-to-noise levels and imaging penetra- 
tion into tissue. 

Clearly, the presented approach requires the use of two or more 
laser sources, depending on the number of acquired wavelengths. 
This comes at the expense of a higher implementation cost and more 
design complexity of the data acquisition system and synchroniza- 
tion electronics. Thus, ideally, alternative methods shall be sought for 
generation of multiple pulses with microsecond delays. While gen- 
eration of low energy pulses in the nanojoule and microjoule range 
with microsecond level delays is indeed possible using pulsed laser 
diodes 28 ' 29 , a different technology will be necessary to manufacture 
affordable devices for deep-tissue multispectral imaging with pulse 
energies in the tens of millijoules range. 

In conclusion, it is anticipated that the presented functional three- 
dimensional optoacoustic imaging approach free of motion artifacts 
will further prompt introduction of the new groundbreaking capabil- 
ities of optoacoustic imaging into clinical practice. 

Methods 

Experimental setup. An optoacoustic hand-held probe based on matrix detection 
array was used to acquire three-dimensional tomographic data in real time 26 . The 
array consists of 256 piezoelectric elements with a central frequency of 4 MHz and a 
bandwidth of ca. 100% ( — 6 dB), which are densely distributed over a spherical 
surface with a radius of 4 cm covering an angle of approximately 90° (solid angle of 
2n[l — cos(tt/4)]). The three-dimensional imaging probe supports nearly isotropic 
spatial resolution of 200 um in a volume of ca. 1 cm 3 surrounding the center of the 
sphere 26 . The output arm of the fiber bundle is inserted into a cylindrical cavity with a 
diameter of 8 mm located in the center of the array to provide optical illumination 
along its central axis. Acoustic coupling to the imaged part is ensured with a water- 
containing transparent polyethylene membrane. Acoustic gel was further applied for 
optimal acoustic coupling between the skin of the imaged subject and the membrane. 
The pressure signals at the 256 locations of the array elements are simultaneously 
collected with a sampling frequency of 40 megasamples per second at 2030 instants, 
immediately following the trigger signal of the first laser pulse. A custom-designed 
data acquisition system (Falkenstein Mikrosysteme GmbH, Taufkirchen, Germany) 
is used for simultaneous digital sampling of the analog data from the 256 channels. 
The acquisition is triggered with the Q- switch output of the first laser, and the 
acquired signals are subsequently transmitted to a personal computer through a 1 
Gbit Ethernet connection. 

Image reconstruction and unmixing. It has been long recognized that extraction of 
quantified information from optoacoustic images is strongly conditioned by accuracy 
of the employed image reconstruction and spectral unmixing schemes 30 . For image 
reconstruction in multispectral optoacoustic tomography, the acquired signals at 
each wavelength are first deconvolved with the electrical impulse response of the 
transducer and band-pass filtered between 0.2 and 7 MHz. An accurate model-based 
reconstruction algorithm 31 is then applied to reconstruct the images for each 
wavelength. When properly reconstructed, the optoacoustic image OA(r,l) y 
representing the pressure field generated by absorption of light energy at a certain 
wavelength X, is proportional to the product of the optical absorption coefficient 
fi a (r,X), the light fluence (j)(r,A) and the dimensionless Griineisen parameter T(r). By 
further considering the optical absorption as superposition of contributions from the 
different absorbing substances present in tissue, OA(r,X) can be expressed as 

OA(r,l) = T(r)(l)(r,X)J2i S iW C M' W 
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where e,-(A) and Q(r) are the molar extinction coefficient and the spatial distribution 
(concentration) of the individual absorbing substances. Assuming that the main 
tissue chromophores in the imaged wavelength range are oxygenated (Hb0 2 ) and 
deoxygenated (HbR) hemoglobin and the main wavelength dependence arises from 
the extinction spectra of hemoglobin, Eq. (1) can be approximated as 



OA(r^T{r)^{r)[8 HbR WC HbR (r)^s Hb o2WC H b02{r)l 



(2) 



where 0(r) is the mean light fluence distribution for the wavelength range employed. 
In order to estimate the concentrations of the two chromophoric components, a 
spectral unmixing procedure based on a per-pixel least square fitting of the 
optoacoustic images at different wavelengths to a linear combination of molar 
extinction coefficients of oxygenated and deoxygenated hemoglobin is performed 32 . 
Subsequently, the spatial distribution of oxygen saturation is estimated via 



S0 2 (r) 



C H boi{r) 



ChwkM + C H bR.{r) 



(3) 



Uncertainty in the oxygen saturation calculation. Since S0 2 (r) is calculated as a 
ratio, it is not directly influenced by quantitative errors in the reconstructed single- 
wavelength images. Nevertheless, the entire procedure relies on approximations 
made in deriving Eq. (2) so that inaccuracies are possibly introduced in cases of strong 
wavelength dependence of the light attenuation, e.g. when considering deep tissue 
imaging scenarios. Other factors causing significant spectral dependence of 
optoacoustic signals, such as wavelength-dependent illuminating beam profile, or 
presence of additional absorbers may all lead to further errors in the estimation of 
C HbR (r) and C Hb02 (r) via Eq. (2). One may estimate the expected inaccuracy of the 
blood oxygen saturation estimation due to wavelength dependence of the light 
attenuation using a recently published procedure 33 . Since in our experiments the 
surface of the imaged object is illuminated with a broad light beam, a first-order 
assessment of the wavelength- dependent light fluence as a function of depth (j)(r,X) 
can be obtained with a one-dimensional light diffusion model considering a uniform 
background 33 . Since the detected optoacoustic signals are normalized to the laser 
energy measured on a per-pulse basis, the light intensity at the air-tissue interface 
(semi- infinite medium) is assumed to be equal for both wavelengths. In calculating 
the wavelength dependent fluence, representative values of the background reduced 
scattering coefficient of fi s '(k) = 18.9(/l/500 nm)" 1 286 and average blood 
concentration of 5.89 g/1 in the forearm were assumed 34 . The background optical 
absorption coefficient /u a for both wavelengths is calculated for this average blood 
concentration and different oxygen saturation levels. Considering a spatially- 
independent Griineisen parameter in soft tissues, the optoacoustically- generated 
signals in each image voxel can be calculated according to Eq. (1) for a given oxygen 
saturation s0 2 in the target voxel and a uniform oxygen saturation s0 2 , b in the 
background tissues. From the generated optoacoustic signals, the approximated 
oxygen saturation s0 2 , a in each voxels can be directly extracted via Eq. (2), which 
assumes the same light fluence distribution for both wavelengths. In this way, the 
expected error in the estimation of the oxygen saturation in the target voxel is given by 



es02 = \s0 2 ,a-s0 2 \ 



(4) 



In vivo human imaging. The imaging experiment in a healthy volunteer was 
performed in full accordance with work safety regulations of Helmholtz Zentrum 
Munchen. The per-pulse energy and average power levels of the lasers conformed to 
standards for the safe use of lasers at all times 27 . To avoid signal saturation due to 
strong optoacoustic responses from the surface, the hair in the imaged region was 
removed prior to the experiments using a shaving lotion. The hand-held probe was 
moved along the skin surface with scanning speeds ranging from 1 to 10 mm/s, as 
estimated in the post-processing step. 
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